stats 153 final project

Importing Data

wholesale_computer = scan("~/desktop/computerwholesale.csv")
wholesale_computer_timeseries = ts(wholesale_computer, frequency = 12, start = c(1997, 1))

Calling the packages needed

library(ggplot2)
library(tseries)

## Warning: package 'tseries' was built under R version 3.3.2

library(TSA)

## Loading required package: leaps

## Warning: package 'leaps' was built under R version 3.3.2

## Loading required package: locfit

## locfit 1.5-9.1    2013-03-22

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-12. For overview type 'help("mgcv-package")'.

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

library(forecast)

## Warning: package 'forecast' was built under R version 3.3.2

## Warning in as.POSIXlt.POSIXct(Sys.time()): unknown timezone 'zone/tz/2018e.
## 1.0/zoneinfo/Asia/Shanghai'

## 
## Attaching package: 'forecast'

## The following object is masked from 'package:nlme':
## 
##     getResponse

library(quantmod)

## Warning: package 'quantmod' was built under R version 3.3.2

## Loading required package: xts

## Warning: package 'xts' was built under R version 3.3.2

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 3.3.2

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: TTR

## Warning: package 'TTR' was built under R version 3.3.2

## Version 0.4-0 included new data defaults. See ?getSymbols.

library(rugarch)

## Warning: package 'rugarch' was built under R version 3.3.2

## Loading required package: parallel

## 
## Attaching package: 'rugarch'

## The following object is masked from 'package:stats':
## 
##     sigma

Initial Plot

plot.ts(wholesale_computer_timeseries)

Decompose data into trend, seasonal and remainder components

fit = stl(wholesale_computer_timeseries, s.window = "periodic")
plot(fit)

Spectral Analysis

raw_spec = spec.pgram(wholesale_computer_timeseries, taper = 0)
plot(raw_spec)

plot(raw_spec, log = "no")

abline(v = 1/7, col = "red")
abline(v = 7/24, col = "red")
abline(v = 9/24, col = "red")
text(x = c(0.2, 0.3, 0.4), labels = c("1/7", "7/24", "9/24"), y = c(1500000, 1700000, 2000000), col = "red")

sp = lm(formula = diff(wholesale_computer_timeseries) ~ cos(2*pi*1/7*seq(from = 1, to = 247)) + sin(2*pi*1/7*seq(from = 1, to = 247)) + cos(2*pi*7/24*seq(from = 1, to = 247)) + sin(2*pi*7/24*seq(from = 1, to = 247)) +cos(2*pi*9/24*seq(from = 1, to = 247)) + sin(2*pi*9/24*seq(from = 1, to = 247)), diff(wholesale_computer_timeseries))

curve(38.7 - 26.334*cos(2*pi*1/7*x) +10.967*sin(2*pi*1/7*x) + 11.275*cos(2*pi*7/24*x) -9.285*sin(2*pi*7/24*x) -2.8 *cos(2*pi*9/24*x) +68.481*sin(2*pi*9/24*x), from = 1, to = 247, ylab = "spectral")
lines(diff(wholesale_computer), col = "blue")

qqnorm(window(rstandard(sp)));qqline(window(rstandard(sp)))

shapiro.test(rstandard(sp))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(sp)
## W = 0.98288, p-value = 0.004535

ARIMA Model

#adf test
adf.test(wholesale_computer_timeseries, alternative="stationary", k=0)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  wholesale_computer_timeseries
## Dickey-Fuller = -3.294, Lag order = 0, p-value = 0.07274
## alternative hypothesis: stationary

adf.test(diff(wholesale_computer_timeseries), alternative = "stationary", k=0)

## Warning in adf.test(diff(wholesale_computer_timeseries), alternative =
## "stationary", : p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(wholesale_computer_timeseries)
## Dickey-Fuller = -19.192, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary

acf(diff(wholesale_computer_timeseries))

pacf(diff(wholesale_computer_timeseries))

eacf(diff(wholesale_computer_timeseries))

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x o o x o o o x x  o  o  o 
## 1 x x x o o x o o o o x  x  o  o 
## 2 x x o o o o o o o o o  x  o  o 
## 3 o x x o o o o o o o o  x  o  o 
## 4 o o x o o o o o o o o  x  o  o 
## 5 x x x o o o o o o o o  o  o  o 
## 6 x o x x o o o o o o x  o  o  o 
## 7 x o o x o o o o o o o  x  o  o

mean(diff(wholesale_computer_timeseries))

## [1] 38.74494

#plot to test for stationary
plot.ts(wholesale_computer_timeseries)

plot.ts(diff(wholesale_computer_timeseries))

#ARIMA model
arima(diff(wholesale_computer_timeseries), order = c(2, 1, 2), include.mean = FALSE)

## 
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(2, 1, 2), include.mean = FALSE)
## 
## Coefficients:
##           ar1      ar2      ma1      ma2
##       -0.6591  -0.2978  -0.5741  -0.4259
## s.e.   0.1466   0.0628   0.1448   0.1444
## 
## sigma^2 estimated as 193392:  log likelihood = -1849.43,  aic = 3706.86

arima(diff(wholesale_computer_timeseries), order = c(0, 1, 3), include.mean = FALSE)

## 
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(0, 1, 3), include.mean = FALSE)
## 
## Coefficients:
##           ma1     ma2     ma3
##       -1.2243  0.1870  0.0373
## s.e.   0.0789  0.1462  0.0837
## 
## sigma^2 estimated as 202130:  log likelihood = -1854.79,  aic = 3715.59

arima(diff(wholesale_computer_timeseries), order = c(1, 1, 3), include.mean = FALSE)

## 
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(1, 1, 3), include.mean = FALSE)
## 
## Coefficients:
##           ar1      ma1      ma2     ma3
##       -0.5310  -0.6605  -0.5587  0.2192
## s.e.   0.1824   0.1756   0.1949  0.0564
## 
## sigma^2 estimated as 198717:  log likelihood = -1852.72,  aic = 3713.44

arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)

## 
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)
## 
## Coefficients:
##           ar1      ar2      ma1     ma2      ma3
##       -0.6641  -0.7899  -0.5453  0.1587  -0.6134
## s.e.   0.0999   0.0912   0.1240  0.1699   0.1274
## 
## sigma^2 estimated as 186715:  log likelihood = -1845.06,  aic = 3700.12

arima(diff(wholesale_computer_timeseries), order = c(3, 1, 3), include.mean = FALSE)

## 
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(3, 1, 3), include.mean = FALSE)
## 
## Coefficients:
##           ar1      ar2     ar3      ma1     ma2      ma3
##       -0.6574  -0.7851  0.0039  -0.5506  0.1621  -0.6115
## s.e.   0.1899   0.1773  0.1021   0.1795  0.1735   0.1584
## 
## sigma^2 estimated as 186717:  log likelihood = -1845.06,  aic = 3702.12

#ARIMA(2, 1, 3) gives the minimum AIC value

tsdiag(arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE))

qqnorm(residuals(arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)));qqline(residuals(arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)))

ARIMA Forecast

fit_arima = stats::arima(wholesale_computer_timeseries, order = c(2, 1, 3))
plot(forecast(fit_arima, h = 30), pch = 19, ylab = "computer sales")

hold = window(ts(wholesale_computer_timeseries), start = 224)
fit_no_holdout = stats::arima(ts(wholesale_computer_timeseries[-c(224:248)]), order = c(2, 1, 2))
plot(forecast(fit_no_holdout, h = 24))
lines(ts(wholesale_computer_timeseries))

Seasonal ARIMA

plot(window(wholesale_computer_timeseries), start = c(1997,1), ylab = "computer sales")

## Warning in plot.window(xlim, ylim, log, ...): "start" is not a graphical
## parameter

## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "start" is
## not a graphical parameter

## Warning in axis(1, ...): "start" is not a graphical parameter

## Warning in axis(2, ...): "start" is not a graphical parameter

## Warning in box(...): "start" is not a graphical parameter

Month = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")
points(window(wholesale_computer_timeseries, start = c(1997, 1)), pch = Month)

auto.arima(wholesale_computer_timeseries)

## Series: wholesale_computer_timeseries 
## ARIMA(3,1,2)(1,0,2)[12] with drift 
## 
## Coefficients:

## Warning in sqrt(diag(x$var.coef)): NaNs produced

##           ar1     ar2    ar3     ma1     ma2     sar1    sma1    sma2
##       -0.7351  -0.811  4e-04  0.4938  0.6274  -0.4239  0.1498  -0.419
## s.e.      NaN     NaN    NaN     NaN     NaN      NaN     NaN     NaN
##         drift
##       35.7920
## s.e.  11.6506
## 
## sigma^2 estimated as 167516:  log likelihood=-1833.77
## AIC=3687.54   AICc=3688.47   BIC=3722.63

#An ARIMA(3,1,2)(1,0,2)[12] is given by auto.arima

plot(diff(diff(wholesale_computer_timeseries),lag=12),xlab='Time',
   ylab='First and Seasonal Difference of computer sales')

acf(as.vector(diff(diff(wholesale_computer_timeseries),lag=12)))

#model
m1_wholesale_computer=arima(wholesale_computer_timeseries,order=c(3,1,2),seasonal=list(order=c(1,0,2),
   period=12))
m1_wholesale_computer

## 
## Call:
## arima(x = wholesale_computer_timeseries, order = c(3, 1, 2), seasonal = list(order = c(1, 
##     0, 2), period = 12))
## 
## Coefficients:
##           ar1      ar2     ar3     ma1     ma2     sar1    sma1     sma2
##       -0.6589  -0.7500  0.0446  0.4389  0.6091  -0.5710  0.3275  -0.3995
## s.e.   0.1860   0.1644  0.1006  0.1749  0.1264   0.2037  0.2004   0.0726
## 
## sigma^2 estimated as 166222:  log likelihood = -1837.29,  aic = 3690.57

#Residuals
plot(window(rstandard(m1_wholesale_computer)), start=c(1997,1), ylab='Standardized Residuals', type='o')

## Warning in plot.window(xlim, ylim, log, ...): "start" is not a graphical
## parameter

## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "start" is
## not a graphical parameter

## Warning in axis(1, ...): "start" is not a graphical parameter

## Warning in axis(2, ...): "start" is not a graphical parameter

## Warning in box(...): "start" is not a graphical parameter

abline(h=0)

acf(as.vector(window(rstandard(m1_wholesale_computer),start=c(1997,1))))

qqnorm(window(rstandard(m1_wholesale_computer),start=c(1997,1)))
qqline(window(rstandard(m1_wholesale_computer),start=c(1997,1)))

shapiro.test(rstandard(m1_wholesale_computer))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(m1_wholesale_computer)
## W = 0.98647, p-value = 0.0192

Seasonal ARIMA Forecasting

hold = window(ts(wholesale_computer_timeseries), start = 212)
fit_no_holdout = stats::arima(ts(wholesale_computer_timeseries[-c(212:248)]), order = c(3, 1, 2), seasonal = list(order = c(1, 0, 2), period = 12))
plot(forecast(fit_no_holdout, h = 36))
lines(ts(wholesale_computer_timeseries))

plot(m1_wholesale_computer,n1=c(1997,1),n.ahead=24,xlab='Year',type = 'o',
   ylab='Computer Sales')

fit_seasonal_ARIMA = forecast(stats::arima(wholesale_computer_timeseries,order=c(3,1,2),seasonal=list(order=c(1,0,2),period=12)))
fit_seasonal_ARIMA$mean

##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 2017                                                               
## 2018 20759.58 20853.65 20535.43 20532.17 20321.25 20333.16 20262.69
## 2019 20075.59 20002.97 20093.43 20015.98 20013.27 20064.48 19854.82
##           Aug      Sep      Oct      Nov      Dec
## 2017          20342.95 20485.66 20868.22 20495.46
## 2018 20249.26 20151.69 20213.96 20142.66 20143.02
## 2019 19583.76

fit_seasonal_ARIMA$lower

##               80%      95%
## Sep 2017 19820.45 19543.86
## Oct 2017 19823.03 19472.26
## Nov 2017 20089.20 19676.81
## Dec 2017 19560.05 19064.87
## Jan 2018 19725.94 19178.76
## Feb 2018 19741.60 19152.91
## Mar 2018 19318.99 18675.04
## Apr 2018 19232.66 18544.74
## May 2018 18956.36 18233.84
## Jun 2018 18888.26 18123.38
## Jul 2018 18744.81 17941.29
## Aug 2018 18672.56 17837.91
## Sep 2018 18539.98 17686.80
## Oct 2018 18560.87 17685.77
## Nov 2018 18455.82 17562.86
## Dec 2018 18424.69 17515.06
## Jan 2019 18320.29 17391.10
## Feb 2019 18214.53 17267.78
## Mar 2019 18274.89 17312.21
## Apr 2019 18163.82 17183.34
## May 2019 18129.04 17131.59
## Jun 2019 18151.11 17138.23
## Jul 2019 17910.36 16881.02
## Aug 2019 17608.53 16562.91

fit_seasonal_ARIMA$upper

##               80%      95%
## Sep 2017 20865.45 21142.04
## Oct 2017 21148.28 21499.06
## Nov 2017 21647.25 22059.64
## Dec 2017 21430.87 21926.05
## Jan 2018 21793.23 22340.40
## Feb 2018 21965.70 22554.39
## Mar 2018 21751.87 22395.81
## Apr 2018 21831.69 22519.61
## May 2018 21686.13 22408.65
## Jun 2018 21778.06 22542.95
## Jul 2018 21780.57 22584.09
## Aug 2018 21825.96 22660.61
## Sep 2018 21763.40 22616.58
## Oct 2018 21867.06 22742.15
## Nov 2018 21829.50 22722.46
## Dec 2018 21861.36 22770.99
## Jan 2019 21830.88 22760.07
## Feb 2019 21791.41 22738.16
## Mar 2019 21911.98 22874.66
## Apr 2019 21868.14 22848.61
## May 2019 21897.51 22894.96
## Jun 2019 21977.84 22990.72
## Jul 2019 21799.28 22828.62
## Aug 2019 21558.98 22604.60

plot(fit_seasonal_ARIMA$residuals)
abline(h = 0)

ARCH/GARCH

diff_wholesales_computer_timeseries = diff(wholesale_computer_timeseries)
plot(diff_wholesales_computer_timeseries);abline(h = 0)

acf(diff_wholesales_computer_timeseries)

pacf(diff_wholesales_computer_timeseries)

acf(abs(diff_wholesales_computer_timeseries))

pacf(abs(diff_wholesales_computer_timeseries))

acf(diff_wholesales_computer_timeseries^2)

pacf(diff_wholesales_computer_timeseries^2)

McLeod.Li.test(y = diff_wholesales_computer_timeseries)

garch_wholesale<-ugarchspec(variance.model = list(model="sGARCH", garchOrder=c(1,1), mean.model = list(armaOrder=c(0,0))), distribution.model = "std")

## Warning: unidentified option(s) in variance.model:
##  mean.model

wholesale_garch<-ugarchfit(spec=garch_wholesale, data=diff_wholesales_computer_timeseries)

wholesale_predict<-ugarchboot(wholesale_garch,n.ahead=10, method=c("Partial","Full")[1])
plot(wholesale_predict,which=2)

ga = garch(diff_wholesales_computer_timeseries, order = c(1,1))

## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     1.929208e+05     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1  1.633e+03
##      1    4  1.633e+03  7.55e-05  2.44e-04  1.1e-07  2.1e+02  4.3e-02  2.57e-02
##      2    5  1.633e+03  2.88e-05  7.54e-05  1.1e-07  2.0e+00  4.3e-02  2.81e-03
##      3    7  1.633e+03  5.97e-06  6.38e-06  8.9e-09  4.5e+00  4.3e-03  6.24e-05
##      4    8  1.633e+03  3.45e-06  4.06e-06  2.2e-08  1.8e+00  8.7e-03  6.86e-06
##      5   10  1.633e+03  1.51e-07  2.77e-07  5.1e-09  5.6e+00  2.2e-03  3.16e-06
##      6   11  1.633e+03  2.55e-08  1.41e-07  5.7e-09  1.3e+00  2.2e-03  1.87e-07
##      7   12  1.633e+03  3.19e-08  3.19e-08  2.1e-09  0.0e+00  8.3e-04  3.19e-08
## 
##  ***** X-CONVERGENCE *****
## 
##  FUNCTION     1.632650e+03   RELDX        2.139e-09
##  FUNC. EVALS      12         GRAD. EVALS       8
##  PRELDF       3.193e-08      NPRELDF      3.193e-08
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    1.929208e+05     1.000e+00     6.439e-06
##      2    9.873919e-02     1.000e+00     8.839e-04
##      3    1.944549e-02     1.000e+00    -1.439e-04

## Warning in garch(diff_wholesales_computer_timeseries, order = c(1, 1)):
## singular information

qqnorm(ga$residuals);qqline(ga$residuals)